Cluster-variation approximation for a network-forming lattice-fluid model 
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We consider a S-dimensional lattice model of a network-forming fluid, which has been recently 
investigated by Girardi and coworkers by means of Monte Carlo simulations [J. Chem. Phys. 126, 
064503 (2007)], with the aim of describing water anomalies. We develop an approximate semi- 
analytical calculation, based on a cluster-variation technique, which turns out to reproduce almost 
quantitatively different thermodynamic properties and phase transitions determined by the Monte 
Carlo method. Nevertheless, our calculation points out the existence of two different phases char- 
acterized by long-range orientational order, and of critical transitions between them and to a high- 
temperature orientationally-disordered phase. Also, the existence of such critical lines allows us to 
explain certain "kinks" in the isotherms and isobars determined by the Monte Carlo analysis. The 
picture of the phase diagram becomes much more complex and richer, though unfortunately less 
suitable to describe real water. 

PACS numbers: 05.50.-|-q 61.20.Gy 65.20.-w 



I. INTRODUCTION 



Water attracts great interest from physicists, due to its 
different anomalous properties [1, 2, 3]. Just to mention a 
few of them, it is well known that, at ordinary pressures, 
the solid phase (ice) is less dense than the correspond- 
ing liquid phase, and the latter displays a temperature 
of maximum density, just above the freezing transition. 
Furthermore, the heat capacity of liquid water is unusu- 
ally large, whereas both isothermal compressibility and 
isobaric heat capacity display a minimum as a function 
of temperature. Although a full prediction of such and 
other anomalies from first principles has not been given 
yet, it is widely believed that it should be related to the 
ability of water molecules to form a network of hydrogen 
bonds. 

Among different possible approaches, a substantial 
body of theoretical investigations concerning so-called 
network-forming fluids has been developed in the frame- 
work of lattice models [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 
14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28]. 
One can use various types of model molecules, charac- 
terized by orientation-dependent interactions, in either 
two [4, 5, 6, 7, 8, 9] or three [10, 11, 12, 13, 14, 15, 16, 
17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28] dimensions. 
A natural choice for water appears to be a 3-dimensional 
model molecule with four bonding arms arranged in a 
tetrahedral symmetry [10, 11, 12, 13, 14, 15, 16, 17, 18, 
19, 20, 21, 22, 23, 24]. Two of them can represent the 
hydrogen (H) atoms, which are positively charged and 
act as donors for the H bond, whereas the other two can 
represent the negatively charged regions ("lone pairs") 
present in the H2O molecule, acting as H bond acceptors. 
As far as the lattice is concerned, the body-centered cu- 



bic (bcc) lattice is suitable for the tetrahedral molecule, 
as the latter can point its arms towards four out of eight 
nearest neighbors of each given site. The above features 
are common to several models, which differ in the form 
of interactions and allowed configurations. 

In the early model proposed by Bell [10, 11, 12], 
molecules can point their arms only to nearest neigh- 
bor sites. An attractive energy is assigned for every pair 
of occupied nearest neighbors, with an extra contribu- 
tion if a hydrogen bond is formed, i.e., if a donor arm 
points to an acceptor arm. Furthermore, a repulsive en- 
ergy is assigned for certain triplets of occupied sites, in 
order to account for the difficulty of forming H bonds by 
closely packed water molecules. Minor variations of the 
Bell model have been proposed by Meijer and cowork- 
ers [13, 14], which replaced the three-body interaction 
with a simple next-nearest-neighbor repulsion, and by 
Lavis and Southern [15], who defined a simplified model 
with no distinction between donors and acceptors, first 
pointing out that such a distinction is scarcely important 
for a qualitative description of water thermodynamics. 
All these studies are mainly focused on the equilibrium 
phase diagram of water, and predict a liquid-vapor co- 
existence and two different low-temperature phases char- 
acterized by long-range orientational order and different 
densities. At zero temperature, the low-density phase is 
an ideal diamond network made up of H bonded water 
molecules, with half the lattice sites left empty, resem- 
bling the structure of ice Ic (cubic ice). At the melting 
temperature, this phase correctly displays a lower density 
than the liquid phase. Conversely, the high-density phase 
at zero temperature is made up of two interpenetrating 
diamond structures, with all sites occupied, resembling 
the structure of ice VII (a high pressure form of ice). 

More recently, a similar model has been proposed by 
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Roberts and Dcbencdetti [18, 19, 20], with the aim of ex- 
ploring the phase diagram of water deep into the super- 
cooled liquid region, and, more generally, the possibility 
of liquid-liquid immiscibility in a network-forming fluid. 
There are two differences with respect to the original Bell 
model. First, water molecules have an additional number 
of configurations, in which they cannot form H-bonds. 
Such a number is an adjustable parameter, related to 
the entropy of breaking H bonds. Second, the energy 
penalty for the triplets occurs only when two molecules 
in a triplet form a H bond. At zero temperature, this 
model still predicts the two ordered H-bond networks of 
the Bell model, but these phases have never been inves- 
tigated at finite temperature, due to the main attention 
devoted by the cited papers to (stable or metastable) 
liquid phases. In particular, these studies have given a 
contribution in the framework of a debate between two 
different conjectures put forth respectively by Speedy and 
Angell [29] and Poole and coworkers [30] to explain ap- 
parent power-law divergences in the response functions 
of supercooled water [29]. In a few words, Speedy and 
Angell postulated a reentrance of the stability limit (spin- 
odal) of liquid water, whereas Poole and coworkers pos- 
tulated the existence of two different metastable liquid 
phases (in analogy with the two different amorphous ices 
observed in experiments), whose coexistence line could 
possibly end in a (metastable) critical point. The model 
by Roberts and Debenedetti turned out to be compati- 
ble with the second critical point conjecture, which has 
also collected several other indirect evidences, both in 
experiments [31] and computer simulations [32]. Two 
authors of the present paper have studied a simplified 
version of the Roberts-Debenedetti model, without the 
donor- acceptor asymmetry [21, 22], pointing out that the 
model is indeed compatible with both the aforementioned 
conjectures (depending on parameter values), and is also 
able to describe the anomalous properties of water as a 
solvent for apolar molecules (hydrophobic effect). 

A different variation of the Bell model has been consid- 
ered by Besseling and Lyklema [16, 17], who have taken 
into account only nearest-neighbor interactions, namely, 
an attractive term for H-bonded molecules and a repul- 
sive term for nonbonded molecules. Even in this case, 
the two different ordered phases described above are sta- 
ble at zero temperature, but they have not been investi- 
gated at finite temperature. The cited studies were in- 
deed focused on the liquid- vapor interface properties [16] 
and on hydrophobic hydration thermodynamics [17], for 
which the authors obtained good agreement with experi- 
ments. These results were obtained in the so-called quasi- 
chemical or Bethe approximation, i.e., a first-order ap- 
proximation which takes into account correlations over 
clusters made up of two nearest-neighbor sites. 

Very recently, Girardi and coworkers [23] have per- 
formed extensive Monte Carlo simulations for the previ- 
ously described model, in the simplified version obtained 
by removing the donor-acceptor distinction. This work 
suggests the existence of two liquid phases of different 



densities, which can coexist in equilibrium, in agreement 
with the conjecture of Poole and coworkers [30]. The 
coexistence line terminates in a critical point, whereas 
the low-density phase exhibits a temperature of maxi- 
mum density (TMD), depending on pressure. A subse- 
quent work [24] suggests a relationship between the den- 
sity anomaly and the anomalous behavior of the diffusion 
coefficient. 

In this paper we analyze the above model by a general- 
ized first-order approximation, based on a 4-site tetrahe- 
dral cluster. Such a work is motivated by the fact that the 
reliability of cluster-variational approximations, though 
widely employed for this kind of models with orientation- 
dependent interactions, has been scarcely tested against 
Monte Carlo simulations [33]. Let us anticipate that 
our calculation reproduces most physical properties ob- 
tained by simulations (namely, isotherms, isobars, bin- 
odals, TMD locus) with remarkable accuracy. Never- 
theless, the picture of the phase diagram turns out to 
be quite different, as the two different condensed phases 
exhibit long-range orientational order, which lead us to 
identify them as a continuous temperature evolution of 
the two zero-temperature network structures introduced 
above. As a consequence, the two critical points are in- 
deed tricritical, and two critical lines appear, indicating 
two different kinds of symmetry breaking, respectively 
from an orientationally-disordered phase to the high- 
density ordered phase and from the latter to the low- 
density one. In order to check the existence of critical 
lines and symmetry breaking we have performed some 
new Monte Carlo simulations in specific regions of the 
phase diagram. 

The paper is organized as follows. In Section II we 
describe the model hamiltonian and analyze the ground- 
state. In Sec. Ill we introduce the cluster-variational 
technique employed for the calculation. Sec. IV describes 
our results, discussing them in comparison with those 
obtained from numerical simulations in Ref. [23]. Sec. V 
reports on the new Monte Carlo results, while Sec. VI is 
devoted to some concluding remarks. 



II. MODEL AND GROUND STATE 

As mentioned in the Introduction, the model is defined 
on a body-centered cubic lattice (see Fig. 1). Each site 
can be empty or occupied by a molecule. The model 
molecule possesses the tetrahedral structure introduced 
above, with 4 equivalent bonding arms, which can point 
towards 4 out of 8 nearest neighbors of a given site. A 
hydrogen bond is formed, yielding an attractive energy 
— ?7 < 0, whenever two nearest-neighbor molecules point 
an arm to each other, with no distinction between donor 
and acceptor (see Fig. 2). Moreover, a repulsive interac- 
tion e > is assigned to any pair of nearest-neighbor sites 
occupied by water molecules. This choice of the interac- 
tion parameters implies an energetic penalty for nearest- 
neighbor molecules not forming H bonds. Finally, as we 
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FIG. 1: Two conventional (cubic) cells of the body centered 
cubic (bcc) lattice. A,B,C,D denote four interpenetrating 
face-centered cubic (fee) sublattices. 




FIG. 2: Two model molecules forming a H bond. The lower 
molecule is in the i = 1 configuration, the upper one is in the 
i = 2 configuration (see the text). 



find it convenient to work in the grand-canonical ensem- 
ble, a chemical potential contribution — ^ is taken into 
account for every occupied site. 

In principle, the hamiltonian of the system can be writ- 
ten as a sum of coupling terms between nearest-neighbor 
sites. Nevertheless, with a view to subsequent analyti- 
cal development, let us write the hamiltonian as a sum 
over irregular tetrahedral clusters, whose vertices lie on 
4 different face-centered cubic sublattices [see Figs. 1 
and 3(a)]. Let us note that there are 24 such tetrahc- 
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FIG. 3: (a) Basic cluster (irregular tetrahedron): A,B,C,D 
denote sites in the 4 corresponding sublattices. AB, BC, 
CD, and DA are nearest-neighbor pairs; AC and BD are 
next-nearest-neighbor pairs, (b) Husimi tree structure corre- 
sponding to the generalized first-order approximation on the 
tetrahedron. 

dra sharing a given site, but it is sufficient to choose only 
a subset of 4, in order to cover all the nearest-neighbor 
pairs [Fig. 3(b)]. There are indeed 6 different possibihties 
to choose a proper subset, but the hamiltonian turns out 
to be independent of this choice. We thus finally write 

W= ^ T-li^ipi^is, (1) 

(a,/3,7,<5) 

where the elementary contribution Tiijki will be denoted 
as tetrahedron hamiltonian. The subscripts ia,ip,i-y,is 
denote site configurations for the 4 vertices a,(3,"f,S of 
each tetrahedron. Possible site configurations are as fol- 
lows: i — denotes a vacancy (empty site), whereas 
J = 1,2 denote a molecule in its 2 possible orientations 
(see Fig. 2). Due to the presence of only nearest-neighbor 
interactions, and assuming that {j,k), {k,l), and 

just refer to nearest-neighbor configurations, the 
tetrahedron hamiltonian can be written as 

Ti-ijki = Ti-ij + Hjk + Ti-ki + Ti-H, (2) 

where 

Hij = eriirij — rjhij — jirii/A. (3) 

In the latter equation, ni is an occupation variable, de- 
fined as rii = if i = (empty site), Ui = 1 other- 
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wise (occupied site), whereas hij is a bond variable, de- 
fined as hij = 1 if the pair configuration repre- 
sents a H bond, and hij ~ otherwise. Let us also as- 
sume that i,j, fc, I (in this order) denote configurations of 
sites placed on, say. A, B, C, D sublattices respectively. If 
A,B,C,D are defined as in Fig. 1, we can define hij = 1 
if i = 1 and j = 2, and hij — otherwise. The 1/4 coeffi- 
cient is meant to avoid multiple countings of the chemical 
potential terms. 

Let us now define the tetrahedron configuration prob- 
ability by Pijkl , with the same convention about the sub- 
script order, and assume that the probability distribution 
is equal for every tetrahedron. Such an assumption seems 
to be reasonable, as it allows one to characterize the two 
ordered structures present in the ground state. In both 
cases, the distribution is a delta function concentrated on 
a fixed tetrahedron configuration k, I) (i.e., pijki ~ 1 
for that configuration, and otherwise), so that there 
is no residual entropy. As far as the low-density (LD) 
structure is concerned, it indeed turns out to be four- 
fold degenerate, since it can be represented by four alter- 
native tetrahedron configurations, namely, {i^j,k,l) = 
(1, 2, 0, 0), (0, 1, 2, 0), (0, 0, 1, 2), (2, 0, 0, 1) (obtained from 
one another by a circular permutation). In each struc- 
ture, a pair of sublattices (respectively, AB, BC, CD, 
and DA) is occupied by fully bonded molecules (form- 
ing a diamond structure), while the other two sublat- 
tices are empty. Conversely, the high-density (HD) struc- 
ture turns out to be two-fold degenerate, being repre- 
sented by the two alternative configurations (i,j,k,l) = 
(f,2,I,2),(2,l,2,l). Both these configurations have all 
lattice sites occupied, but they differ in the pairs of fully 
bonded sublattices, that is, AB and CD in the former 
case and BC and DA in the latter. 

The grand-canonical energy per site in the thermody- 
namic limit can be written as 

2 

w= ^ PijkiTi-ijki, (4) 
ij,k,l=0 

which takes into account that there is one tetrahedron 
per site. By grand-canonical energy, we mean 

w = u — fip, (5) 

where u is the internal energy per site and p is the density, 
i.e., the average occupation probability 

2 

P^jkl ■ (6) 

We arc now in a position to evaluate the energy of the 
two ground-state structures described above. For the LD 
phase, we have to replace 

Pijki = <5i,i(5j,2(5fc,o'5i,o (7) 

into Eq. (4), obtaining 

WLD =e-ri- fi/2, (8) 



while Eq. (6) correctly yields the density pld = 1/2. 
Moreover, for the HD phase, we consider 

Pijkl = Si^i6j^2Sk,iSi^2 (9) 

which; replaced into Eqs. (4) and (6), yields respectively 

wuD = 4:6 - 2r] - p, (10) 

and phd — 1- Finally, wc have also to consider a gas 
(G) phase, in which all but a finite number of sites are 
empty, so that both density and energy are in the ther- 
modynamic limit. It is easy to show that the G phase is 
thermodynamically favored (0 < wld and < whd) for 
p < ^G-LD, where 

MG-LD = 2e - 27/, (II) 

the LD phase is favored (wld < and wld < whd) for 
MG-LD < ^J. < MLD-HD, where 

Mld-hd = 6e - 2r], (12) 

and the HD phase is favored (whd < and whd < u'ld) 
for ^ > ^LD-HD- Let us note that the LD phase has 
always a stability region, as /xg-ld < Mld-hDj because of 
the repulsive term e > 0. 



III. CLUSTER- VARIATION APPROXIMATION 

We have performed the finite temperature analysis 
by means of a cluster-variational approximation, previ- 
ously applied by two of the authors for investigating a 
similar water-like lattice model [21, 22]. The cluster- 
variation method is an improved mean-field theory, which 
in principle can take into account correlations at arbitrar- 
ily large, though finite, distances. In Kikuchi's original 
work [34], an approximate entropy expression was ob- 
tained by heuristic arguments, while, in more recent and 
rigorous formulations [35] , the approximation is shown to 
be equivalent to a truncation of a cluster cumulant ex- 
pansion of the entropy. The approximation is expected 
to work, because of a rapid decreasing of the cumulant 
magnitude, upon increasing the cluster size, namely when 
the latter becomes larger than the correlation length of 
the system [36]. A particular approximation is defined 
by the largest clusters left in the truncated expansion, 
usually denoted as basic clusters. One obtains a free en- 
ergy functional in the cluster probability distributions, 
to be minimized, according to the variational principle 
of statistical mechanics. 

For our model we have chosen, as basic clusters, the 
irregular tetrahedra which we have previously used to 
express the model hamiltonian Eq. (I). Indeed, this 
choice turns out to coincide with the (generalized) first- 
order approximation (on the tetrahedron cluster), which 
is also equivalent to the exact calculation for a Husimi 
lattice [37], whose (tetrahedral) building blocks arc just 
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arranged as in Fig. 3(b). The main advantage of this ap- 
proximation is its simpUcity, as the only clusters retained 
in the expansion are basic clusters and single sites (it is 
sometimes denoted as cluster-site approximation [38]), 
together with a substantial improvement over the ordi- 
nary mean-field theory, which has been recognized for 
different models, even with orientation dependent inter- 
actions [9]. In particular, let us note that the internal 
energy is treated exactly, since the range of interactions 
is smaller than the basic cluster size. The grand canonical 
free energy per site u) = w — Ts (T being the tempera- 
ture, expressed in energy units, and s the entropy per 
site, in natural units) can be written as 



OJ 

T 



2 

i,j,k,l=0 



T 



In 



Pijkl 



[pfpfp'^pF] 



(13) 

where pf is the probability of the i configuration for a 
site on the X sublattice {X = A, B, C, D). These proba- 
bilities can be obtained as marginals of the tetrahedron 
probability distribution as 



Pijkl, 



pf ^ 

j,k,i=a 

2 





k,lA=0 
2 



(14) 



Accordingly, the variational free energy in Eq. (13) turns 
out to be a function of only the tetrahedron distribution 
Pijkl ■ Let us note that the free energy expression contains 
two different logarithmic terms, respectively correspond- 
ing to the cluster and site entropies. The former takes 
into account correlations among four configuration vari- 
ables on a tetrahedron, while the latter can be viewed as a 
correction term ensuring that, if the tetrahedron distribu- 
tion factorizes into a product of single site probabilities, 
the mean field (Bragg- Williams) entropy approximation 
is recovered. By the way, such a heuristic argument was 
first used by Guggenheim [39] to derive an equivalent ex- 
pression for the case of a two-site cluster, also equivalent 
to the well-known Bethe approximation [40]. 

Minimization of lo with respect to the set of tetrahe- 
dron probabilities {pijki\, with the normalization con- 
straint 



i,j,k,l=G 



(15) 



can be performed by the Lagrange multiplier method, 
yielding the equation 

K,M=r^e-«'-'/^(p^fp^pf)-^/\ (16) 

where ^ is related to the Lagrange multiplier, and can be 
obtained imposing the constraint Eq. (15) 

e= E e-^^^^^'^ [pfpfp'iprf'. (17) 

i,j,k,l=Q 



Eq. (16) is in a fixed point form, and can be 
solved numerically by simple iteration (natural iteration 
method [41]). For the cluster-site approximation, the nu- 
merical procedure can be proved to reduce the free energy 
at each iteration [37, 41], and therefore to converge to lo- 
cal minima. From the solution of Eq. (16) one obtains 
a tetrahedron distribution pijku whence one can com- 
pute the thermal average of every observable; density by 
Eq. (6), internal energy by Eqs. (4) and (6), and free 
energy by Eq. (13). The latter can be also related to the 
normalization constant as [37] 



-Tln^. 



(18) 



Finally, assuming the volume per site equal to unity, pres- 
sure can be determined, in energy units, as P = —oj. 
In the presence of multiple solutions, i.e., of competing 
phases, the thermodynamically stable one is selected by 
the lowest free energy (highest pressure) value. At zero 
temperature, we have P = —w, so that we can easily 
determine the transition pressures as follows. Replacing 
Eq. (11) into (8), we obtain the G-LD transition pressure 



G-LD 



0, 



(19) 



while, replacing Eq. (12) into (10), we obtain the LD-HD 
transition pressure 

Plb-hb = 2e. (20) 

IV. RESULTS AND DISCUSSION 

In this section we present and discuss our results, in- 
cluding a detailed comparison with the Monte Carlo sim- 
ulations of Girardi and coworkers [23]. With this pur- 
pose, we consider the same ratio of interaction parame- 
ters as in Ref. 23, namely, 7]/e = 2. Let us note that the 
ground-state transition pressures computed above are in- 
dependent of this parameter. 

A. Phase diagram 

In Figs. 4 and 5 we report two projections of the phase 
diagram, respectively in the temperature-pressure and 
density-temperature planes. We can observe three differ- 
ent phases, which we denote as G, LD, and HD, as they 
can be respectively identified as continuous evolutions 
of the three ground-state phases discussed above. The 
G phase is homogenous, in that the single-site probabil- 
ity distribution is independent of the sublattice, whereas 
the LD and HD phases exhibit a symmetry-breaking, as 
different sublattices behave differently. In the homoge- 
neous phase, molecules can locally form H bonds, but 
a long-range ordered H-bond network does not appear. 
Conversely, the other two phases exhibit long-range or- 
der, and one can identify the same types of symmetry- 
breaking observed in the zero-temperature LD and HD 
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FIG. 4: Pressure {P/^) vs temperature {T/e) phase diagram. 
G, LD, HD denote the corresponding phase regions (see the 
text); TP and TP' are tricritical points; CEP is the critical 
end-point. Solid and dashed lines denote respectively first- 
and second-order transitions; a dash-dotted line denotes the 
TMD locus. Squares and circles display Monte Carlo data 
from Ref. 23 for the first-order transitions and the TMD locus, 
respectively. 




FIG. 5: Temperature {T/e) vs density (p) phase diagram. La- 
bels are defined as in Fig. 4 (double labels denote coexistence 
regions). Solid lines denote boundaries of the coexistence re- 
gions, while dashed lines denote second-order transitions; a 
dash-dotted line denotes the TMD locus. Squares and circles 
display Monte Carlo data from Ref. 23 for the coexistence 
boundaries and the TMD locus, respectively. 



phases, respectively. More precisely, upon decreasing 
temperature, the tetrahedron distributions tend to be 
concentrated on the fixed configurations, corresponding 
respectively to the different (LD or HD) H-bond networks 
described above. For instance, in the LD phase, H bonds 
arc preferentially formed through AB (or BC, CD, DA) 
sublattices, whereas, in the HD phase, through AB and 
CD (or BC and DA). 

Let us have a closer look at the phase diagram topol- 



ogy. In the pressure vs temperature diagram (Fig. 4), we 
observe three different first-order transition lines. The 
first one separates the G and LD phases, whereas the 
other two occur between the LD and HD phases. All 
these lines are mapped onto coexistence regions in the 
temperature vs density diagram (Fig. 5). Both the LD- 
HD first-order transition lines terminates in tricritical 
points (denoted as TP and TP'), which turn out to be 
connected by a second-order transition line (i.e., a line of 
critical points), which encloses the LD phase. Another 
critical line separates the G phase from the HD phase, 
which correctly prevents any continuous transformation 
between these two phases, the latter having a lower sym- 
metry. This critical line terminates in a critical end-point 
(CEP). Let us note, by the way, that, in the LD-HD co- 
existence region between CEP and TP, the LD phase 
has in fact a higher density than HD, so that the two 
phases are identified on the basis of their different sym- 
metries. Let us also note that the LD phase exhibits 
a density anomaly, namely, a temperature of maximum 
density (TMD), at constant pressure. 

Figs. 4 and 5 also display some data point obtained 
by the Monte Carlo analysis of Ref. [23]. We find a 
remarkably good agreement both for first-order transi- 
tions (coexistence regions in Fig. 5 and transition lines 
in Fig. 4), and for the TMD locus in the LD phase. De- 
viations are slightly larger only in the vicinity of critical 
lines, where the correlation length of the system diverges, 
so that we expect a break-down of our approximation 
performances. In spite of such an impressive agreement 
between Monte Carlo results and our approximate calcu- 
lation, which suggests that the tetrahedron approxima- 
tion captures the most relevant correlations present in 
the system, an important difference emerges between the 
two studies. Indeed, Girardi and coworkers [23] recognize 
the two long-range ordered network structures only in the 
ground-state analysis. Such a long-range order seems to 
be lost at finite temperature, so that they denote the 
two denser phases as low density liquid (LDL) and high 
density liquid (HDL). In this scenario, they only observe 
two first-order transitions lines (G-LDL and LDL-HDL), 
terminating in two different critical points. 

Although in principle such a discrepancy might be 
an artifact of our approximate method, different evi- 
dences lead us to conjecture that this is not the case, 
and that in Ref. 23 the authors have simply not ex- 
plored the possibility of a symmetry-breaking. A first 
evidence is of course the striking quantitative agreement 
shown above. Moreover, as previously mentioned, our 
symmetry-broken phases at T > are continuous evolu- 
tions of those predicted at T = 0, whereas in Ref. 23 it is 
not clear how the transitions from the zero-temperature 
ordered phases to the finite temperature liquid phases 
occur. Looking for further evidences, we have analyzed 
isotherms and isobars. 
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B. Isotherms and isobars 

In Fig. 6 wc report isotherms in a pressure vs den- 
sity diagram, together with corresponding Monte Carlo 
data, available from Ref. 23. In particular, we consider 
three different temperatures, showing qualitatively differ- 
ent behaviors. The first one is T/e = 0.8 (Fig. 6a), where 
the isotherm displays two large plateaus, resulting from 
the two different (G-LD and LD-HD) first-order tran- 
sitions discussed above. These plateaus can be clearly 
observed even from the Monte Carlo data, which are 
very well fitted by our curve. The second temperature 
is T/e = 1.2 (Fig. 6b), lying between the two tricritical 
points TP and TP' (see Figs. 4 and 5). In this case, 
upon increasing pressure, the system first undergoes a 
G-LD first-order transition, and then a LD-HD second- 
order transition. The former results in a plateau as well, 
whereas the latter gives rise to a simple "kink" in the 
isotherm. The Monte Carlo data exhibit a similar plateau 
and, noticeably, they also seem to be compatible with the 
possibility of a critical-line crossing, i.e., with the exis- 
tence of a kink in the thermodynamic limit. Finally, we 
consider T/e = 1.6 (Fig. 6b), lying above both tricritical 
points. Here the system directly undergoes a (second- 
order) G-HD transition, so that the isotherm displays 
just a slight kink. Even in this case, the slope change 
observed in the Monte Carlo isotherm seems to be well 
explained by the critical-line crossing. 

Let us note that, looking at Figs. 4 and 5, we indeed 
expect two more qualitatively different regimes, respec- 
tively corresponding to the narrow temperature intervals 
between CEP and TP and between TP and the maximum 
temperature reached along the TP-TP' critical line. In 
the former case, upon increasing pressure, we find a G- 
HD second-order transition, followed by a HD-LD first- 
order transition, and finally by a LD-HD second-order 
transition. In the latter case, we find the same sequence 
of transitions, all second-order. Corresponding isotherms 
are not displayed in Fig. 6. Indeed, we suppose that a 
comparison with Monte Carlo simulations may be very 
difficult in these temperature region, where maximum 
quantitative discrepancies are observed in the phase di- 
agram (see Fig. 5). To be more precise, it is likely that 
similar behaviors might also appear from Monte Carlo 
data, but at a slightly lower temperature than observed 
from the approximate calculation. 

In Fig. 7 we report isobars in a density vs temper- 
ature diagram (more precisely, Fig. 7a displays only re- 
sults from our calculation, whereas in Fig. 7b some Monte 
Carlo results are superimposed). We have considered dif- 
ferent pressure values, below the TP' pressure, pointing 
out four different regimes. The lowest pressure isobars lie 
below the CEP pressure PcEp/e ~ 1-4, so that the sys- 
tem only undergoes a first-order G-LD transition, upon 
decreasing temperature (see Fig. 4). At the lowest pres- 
sures, the density maximum turns out to be extremely 
shallow, but it becomes more pronounced upon increas- 
ing pressure (the TMD locus correctly connects all den- 



sity maxima). At pressure values lying between Pcep 
and the TP pressure Ptp/c ~ 1-8, the system phase be- 
havior becomes more complicated. Upon decreasing tem- 
perature, we first observe a second-order G-HD transition 
followed by a weak first-order HD-LD transition. The 
HD-LD transition also becomes continuous above Ptp- 
Finally at P/e > 2 [i.e., above the ground-state LD-HD 
transition pressure, see Eq. (20)], a further LD-HD first- 
order transition appears at low temperature. A compari- 
son of these results with Monte Carlo isobars drawn from 
Ref. 23 generally confirms the good agreement obtained 
for phase diagrams and isotherms. In particular, let us 
note that the G-HD critical line provides a clear interpre- 
tation of the slope changes observed in the two highest 
pressure isobars in Fig. 7b. 

Looking at Fig. 4, we expect two more regimes, for 
pressure values above TP'. More precisely, between TP' 
and the maximum pressure value reached along the TP- 
TP' critical line, we observe a sequence of three contin- 
uous transitions, namely, G-HD, HD-LD, and LD-HD, 
upon decreasing temperature. For higher pressure, we fi- 
nally observe a direct G-HD continuous transition. Even 
in these cases, corresponding isobars have not been re- 
ported in Fig. 7. 



V. A MONTE CARLO TEST 

As pointed out in the previous section, different ev- 
idences support the existence of long-range ordered 
symmetry-broken phases, and related critical transitions. 
Nevertheless, due to the approximate nature of our re- 
sults, we have found it necessary to perform some specific 
checks by Monte Carlo simulations. We have only inves- 
tigated certain regions around which the transitions are 
expected to appear. 

In order to verify the G-HD transition, we consider a 
constant temperature T/e = 1.7 (see Figs. 5 and 4), and 
work out a grand-canonical simulation for varying chem- 
ical potential. The lattice is made up of T x L x T bcc 
conventional cells (see Fig. 1), with L ranging from 10 up 
to 30, and periodic boundary conditions. The simulation 
is of the Metropolis type, starting from a random lattice 
configuration. We randomly select a lattice site, and fiip 
its configuration according to the Metropolis rule. We 
define the Monte Carlo step as 2L'^ iterations of this pro- 
cess. Measurements are performed every 30-100 Monte 
Carlo steps, depending on the phase (the Monte Carlo 
dynamics is faster in the disordered G phase than in the 
LD and HD phases). We collect up to 10'' measurements, 
waiting a thermalization time of order 10'^ measurements. 

In Fig. 8 wc report the fourth-order cumulant [42] of 
the grand-canonical energy Vw as a function of the chem- 
ical potential, for different values of the lattice size L. We 
can observe that Vw fluctuates around zero, except in a 
narrow interval (which shrinks upon increasing L) where 
Vw exhibits a minimum followed by a maximum. This is a 
typical signature of a second-order transition [43] , which 
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FIG. 6: Isotherms in the pressure (P/e) vs density (p) diagram. Solid and dashed (thick) lines denote respectively boundaries 
of coexistence regions and critical lines (some lines are reported partially, to improve readability). Thin lines of different types 
represent isotherms at T/e = 0.8 (a) and T/e = 1.2, 1.6 (b). Scatters denote Monte Carlo results from Ref. 23. Labels as in 
Fig. 4. 




FIG. 7: Isobars in the density (p) vs temperature (T/e) diagram. Solid and dashed (thick) lines denote respectively boundaries 
of coexistence regions and critical lines. Thin lines of different types represent isobars at various pressures. In (a), a thick 
dash-dotted line denotes the TMD locus. In (b), scatters denote Monte Carlo results from Ref. 23. Labels as in Fig. 5. 



in our case we expect to be the G-HD transition. Actu- 
ally, in both phases far from the transition, the energy 
distribution can be described by a single peak, which can 
be approximated by a gaussian, tending to a Dirac delta 
in the L —f oo limit. In this case Vw should vanish inde- 
pendently of L. Upon approaching the critical line, the 
symmetric gaussian peak description is no longer valid, 
so that Vw ^ [43] . Though we do not report the details 
here, we have worked out an equivalent analysis for the 
LD-HD critical line, obtaining similar results. 

The results displayed in Fig. 8 seem to indicate a tran- 
sition around 1.3 < fj,/e < 1.4, where the simulation 
predicts a density 0.59 ^ p < 0.61. The latter val- 
ues turn out to be slightly higher than the transition 
density predicted by the cluster-variation approximation 
p sa 0.56, at the corresponding transition chemical poten- 
tial /.t/e « 1.00. We then expect that the actual critical 



line is shifted towards higher density (i.e., higher pres- 
sure) with respect to the approximate line shown in Figs. 
4 and 5. This result is in agreement with the known be- 
havior of cluster-variational approximations, which, due 
to their mean-field nature, usually overestimate the or- 
dered phase region. 

In order to better characterize the transition, we define 
the following order parameter 

(21) 

where is the probability of the i — 1 configuration 
(see Fig. 2) on the X and Y sublattices. A totally equiv- 
alent order parameter 02 could be obtained considering 
the i — 2 configuration. Let us note that, as previously 
mentioned, the HD phase is two-fold degenerate. There- 
fore, we expect that our simulated system can go into 
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bilities as 



FIG. 8: Fourth-order cumulant of the grand canonical en- 
ergy (Vw) vs chemical potential (yu/e) at constant tempera- 
ture T/e — 1.7, for different values of the lattice size L. Error 
bars, not reported for a better readability, are of order 10~^ 
in the worst case. 



X 



(22) 



where the one-sublattice probabilities are derived from 
Eqs. (14). It is easy to see that (f>i behaves like a good 
order parameter for the G-HD transition. Indeed, for low 
/i values, where we expect to observe the disordered G 
phase, (/ii tends to zero upon increasing L. Conversely, 
for higher fi values, where we expect to observe the or- 
dered HD phase, (jji remains nonzero, with a smaller ef- 
fect of L. In the intermediate region, the slope of the 
curves increases upon increasing L, revealing the phase 
transition. Let us note, by the way, that for L = 10 (the 
size used in Ref. 23) finite size effects are definitely not 
negligible in the vicinity of the transition, where the cor- 
relation length of the system is expected to diverge. Ac- 
cordingly, also the cluster-variational result fails in this 
region, although it progressively improves upon moving 
away from the transition. 
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FIG. 9: G-HD order parameter (l^il) vs chemical potential 
(/i/e) at constant temperature T/e = 1.7. The solid line de- 
notes the cluster-variation result, while scatters denote Monte 
Carlo results for different values of the lattice size L (thin lines 
are an eyeguide). Error bars are of order lO"**. 



one of the two possible states, depending on initial con- 
ditions. In one case, preferential H bonding occurs cither 
through AB and CD sublattices, with preferred i = 1 
configurations on A and C. In the other case, preferen- 
tial bonding occurs through BC and DA, with preferred 
i = 1 configurations on B and D. In the two cases we 
respectively obtain positive and negative values of the 
order parameter. 

In Fig. 9 we report the absolute value of 4>i , computed 
both by simulations (for different lattice sizes L) and 
by the cluster-variation approximation (in the thermo- 
dynamic limit), as a function of the chemical potential. 
In the latter case, we compute the two-sublatticc proba- 



VI. CONCLUSIONS 

In this paper, we have considered a lattice model with 
orientation-dependent interactions, meant to describe a 
so-called network-forming fluid. As mentioned in the In- 
troduction, this model is an instance of quite a large 
class of similar models (defined on the body-centered 
cubic lattice, with tetrahedral model molecules), origi- 
nally conceived to describe water and investigate its var- 
ious anomalies. Several studies have been carried out 
using a variety of approximate techniques, whose relia- 
bility has not always been checked. Motivated by a re- 
cent work by Girardi and coworkers [23] , performing ex- 
tensive Monte Carlo simulations for one of such models, 
we have worked out a corresponding analysis by means 
of a cluster-variational technique, with the aim of com- 
paring the results and evaluating the reliability of the 
approximate method. The latter, based on a tetrahedral 
cluster, had been already employed by two authors of 
the present paper for investigating a different water-like 
lattice model [21], and actually turns out to be a general- 
ization of the original first-order approximation worked 
out by Bell [10]. The generalization consists in taking 
into account the possibility of a symmetry-breaking at 
the level of four face-centered cubic sublattices of the 
original body-centered cubic lattice. Such a symmetry 
breaking is indeed known to occur in the model ground- 
state. 

As already pointed out in the text, we obtain a re- 
markably good agreement with Monte Carlo results for 
most thermodynamic quantities. Nevertheless, our re- 
sults clearly suggest that the the two denser phases, iden- 
tified as liquid in Ref. 23, are in fact endowed with a long 
range ordered structure, closely resembling the hydrogen 
bond networks observed at zero temperature. Accord- 
ingly, the emerging picture of the phase diagram turns 
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out to be much more complex and richer. In fact, due to 
symmetry-breaking, the two critical points observed in 
Ref. 23 turns out to be tricritical, and two new critical 
lines emerge, separating the two symmetry-broken phases 
from each other and from the disordered phase. Such 
critical lines arc hardly detectable by observing macro- 
scopically averaged quantities (such as the density), since 
they manifest themselves only as kinks (not plateaus) in 
isotherms and isobars, which are further smoothed out 
by finite size effects. We have indeed performed a Monte 
Carlo simulation, in order to check the qualitative cor- 
rectness of our results, around the critical lines deter- 
mined by the approximate method. The density exhibits 
in fact a smooth behavior, but the analysis of a Binder 
cumulant [43], and of a suitable order parameter, shows 



clear evidence of an order-disorder critical transition. 

As previously mentioned, these results suggest that the 
tetrahedral cluster approximation is able to capture the 
most relevant correlations present in the system. On 
the other hand, unfortunately, the observed phase dia- 
gram seems to make the model less relevant for investi- 
gating water. Concerning this issue, the main problem 
is of course that the two ordered phases, which indeed 
may well represent two different ice forms [11], undergo 
second-order transitions, which have never been observed 
experimentally. It is likely that further variations of the 
model might improve the similarity of its phase diagram 
to that of real water, but a discussion of this point is be- 
yond the scope of the present paper and is left for future 
work. 
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